# 0028	Environmental Protection Specialist	Administrative (A)	Inspectors focusing on environmental compliance, permits, enforcement.
# 0819	Environmental Engineer	Professional (P)	Engineering-based inspections, technical assessments.
# 1301	Physical Scientist	Professional (P)	Scientific inspections, environmental sampling, site investigations.
# 1801	General Inspection, Investigation, Enforcement	Administrative (A)	Broader compliance and enforcement inspectors.


# Load necessary libraries
require(tidyverse)  # For data manipulation and visualization
require(sf)

# Load the data
load("00_data/01_data_processed/epa_empl.Rdata")

# Define EPA treatment regions
EPA_TREAT <- c("01", "02", "03", "04", "05", "07")

# Map states to EPA regions
state_to_region <- list(
  `01` = c("CONNECTICUT", "MAINE", "MASSACHUSETTS", "NEW HAMPSHIRE", "RHODE ISLAND", "VERMONT"),
  `02` = c("NEW JERSEY", "NEW YORK", "PUERTO RICO", "U.S. VIRGIN ISLANDS"),
  `03` = c("DELAWARE", "DISTRICT OF COLUMBIA", "MARYLAND", "PENNSYLVANIA", "VIRGINIA", "WEST VIRGINIA"),
  `04` = c("ALABAMA", "FLORIDA", "GEORGIA", "KENTUCKY", "MISSISSIPPI", "NORTH CAROLINA", "SOUTH CAROLINA", "TENNESSEE"),
  `05` = c("ILLINOIS", "INDIANA", "MICHIGAN", "MINNESOTA", "OHIO", "WISCONSIN"),
  `06` = c("ARKANSAS", "LOUISIANA", "NEW MEXICO", "OKLAHOMA", "TEXAS"),
  `07` = c("IOWA", "KANSAS", "MISSOURI", "NEBRASKA"),
  `08` = c("COLORADO", "MONTANA", "NORTH DAKOTA", "SOUTH DAKOTA", "UTAH", "WYOMING"),
  `09` = c("ARIZONA", "CALIFORNIA", "HAWAII", "NEVADA", "AMERICAN SAMOA", "GUAM", "NORTHERN MARIANA ISLANDS"),
  `10` = c("ALASKA", "IDAHO", "OREGON", "WASHINGTON")
)

# Function to find the EPA region for a given state
find_region <- function(state_name) {
  for (region in names(state_to_region)) {
    if (state_name %in% state_to_region[[region]]) {
      return(region)
    }
  }
  return(NA)  # Return NA if the state name does not match any region
}

# Summarize EPA employment data by state and year
EPA_EMPL_COUNTY <- EPA_EMPL_COUNTY %>% 
  filter(occupation %in% c("0028", "0819", "1301", "0801")) %>% 
  group_by(State, year) %>% 
  summarize(
    number = n(),
    median_grade = median(grade_numeric, na.rm = TRUE),
    median_edu = median(education_numeric, na.rm = TRUE),
    mean_grade = mean(grade_numeric, na.rm = TRUE),
    mean_edu = mean(education_numeric, na.rm = TRUE)
  )


# Assign EPA regions to each entry
EPA_EMPL_COUNTY$epa_region <- sapply(EPA_EMPL_COUNTY$State, find_region)

# Indicate if a state is in Phase 1 and if it is in a treatment region
EPA_EMPL_COUNTY$treat_region <- EPA_EMPL_COUNTY$epa_region %in% EPA_TREAT 
EPA_EMPL_COUNTY$time_to_treat <- EPA_EMPL_COUNTY$year - 1995 

# Filter data for treatment years and exclude DC
treat_year <- EPA_EMPL_COUNTY %>% 
  filter(State != "DISTRICT OF COLUMBIA") %>% 
  filter(year >= 1989 & year <= 2001)

# Summarize data by EPA region, year, treatment status, and time to treatment
treat_year_epa <- treat_year %>% 
  group_by(epa_region, year, treat_region, time_to_treat) %>%
  summarize(number = sum(number))


# Plot the percentage change in EPA employees by treatment status over time
treat_year_epa %>% 
  mutate(Treatment = factor(if_else(treat_region == TRUE, 1, 0))) %>% 
  group_by(epa_region) %>% 
  arrange(year) %>% 
  mutate(pct_change = (number / lag(number) - 1) * 100) %>% 
  filter(!is.na(pct_change)) %>% 
  group_by(Treatment, year) %>% 
  summarize(mean_pct_change = mean(pct_change)) %>% 
  ggplot(aes(x = year, y = mean_pct_change, color = Treatment, fill = Treatment)) +
  geom_point() +
  geom_line() +
  theme_bw() +
  scale_x_continuous(breaks = 1990:2001) +
  ylab("Pct. change of EPA employees") +
  xlab("Year") +
  ylim(-10, 20) +
  geom_hline(yintercept = 0, linetype = 1) +
  geom_vline(xintercept = 1994, linetype = 2) +
  geom_vline(xintercept = 1999, linetype = 2) +
  scale_color_manual(values = rev(c("#0A2463", "#FB3640"))) +
  scale_fill_manual(values = rev(c("#0A2463", "#FB3640"))) +
  theme(legend.position = c(0.9, 0.2),
        legend.background = element_rect(fill = "white",
                                         linewidth = 0.2, 
                                         linetype = "solid", 
                                         colour = "black")) +
  coord_cartesian(ylim = c(-20, 20))

# Save the plot to a file
ggsave("03_figures/appendix/figure_a13.pdf", height = 3, width = 8)

rm(list = ls())

